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Abstract. A “dimensional reduction” (DR) method is introduced for analyzing lightning field 
changes (A Es) whereby the number of unknowns in a discrete two-charge model is reduced 
from the standard eight ( x,y,z,Q,x',y',z',Q ') to just four (x,y,z,Q). The four unknowns 
( x,y,z,Q ) are found by performing a numerical minimization of a chi-square function. At each 
step of the minimization, an Overdetermined Fixed Matrix (OFM) method is used to 
immediately retrieve the best “residual source” (x',y',z',Q ') , given the values of ( x,y,z,Q ). In 
this way, all 8 parameters ( x,y,z,Q,x',y',z',Q ') are found, yet a numerical search of only 4 
parameters ( x,y,z,Q ) is required. The DR method has been used to analyze lightning-caused 
A Es derived from multiple ground-based electric field measurements at the NASA Kennedy 
Space Center (KSC) and USAF Eastern Range (ER). The accuracy of the DR method has been 
assessed by comparing retrievals with data provided by the Lightning Detection And Ranging 
(LDAR) system at the KSC-ER, and from least squares error estimation theory, and the method 
is shown to be a useful “stand-alone” charge retrieval tool. Since more than one charge 
distribution describes a finite set of A Es (i.e., solutions are non-unique), and since there can 
exist appreciable differences in the physical characteristics of these solutions, not all DR- 
solutions are physically acceptable. Hence, an alternative and more accurate method of analysis 
is introduced that uses LDAR data to constrain the geometry of the charge solutions, thereby 
removing physically unacceptable retrievals. The charge solutions derived from this method are 
shown to compare well with independent satellite- and ground-based observations of lightning in 
several Florida storms. 
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1. Introduction 


Lightning-caused changes in the surface electric field have been used to determine the 
characteristics (magnitude, polarity, location) of the charges that are effectively deposited by 
lightning [e.g., Wilson, 1916, 1920; Workman and Holzer, 1942; Jacobson and Krider, 1976; 
Krehbiel et ah, 1979; Krehbiel, 1981, 1986; Maier and Krider, 1986; Koshak and Krider, 1989, 
1994; Murphy et ah, 1996; Koshak et ah, 1999]. Conventional charge retrieval algorithms 
minimize a nonlinear, chi-square function to find the optimum parameters, e.g. charge locations 
and values, of an assumed model of the source. Retrieval accuracy depends in part on the model 
and on the error in the lightning-caused electric field changes ( AEs ) . For the ground-based 
network of electric field sensors (field mills) at the NASA Kennedy Space Center (KSC) and the 
USAF Eastern Range (ER), and for the algorithm used to compute the A Es from the electric 
field records (e.g., see Koshak and Krider, 1 989), the estimated error in AE is typically under 
15%. Normally, a simple charge model is used to fit the A Es to within this error, and if this 
approach is successful, it complies with a principle known as Occam’s Razor. The 14 th century 
English philosopher and Franciscan monk William of Occam, has stated " Pluralitas non est 
ponenda sine neccesitate" or "plurality should not be posited without necessity." From this point 
of view, a simple model composed of just one or two spherically symmetric sources would be 
better than a more complicated solution if the simple solution described the observed AEs to 
within 15%. In effect, Occam’s Razor uses only the simplest models to describe the observations 
and ignores all other possibilities. 

In this paper, section 2 briefly reviews the goodness-of-fit method that is customarily 
used to determine the optimum parameters of a charge model, and includes a description of the 
one-charge model (also called the point-charge or “Q-model”), and the two-charge model. 



Section 3a reviews a 4-station analytic method that can be used to find the parameters of the one- 
charge model (Krehbiel et al., 1979), section 3b generalizes this method to 4 or more stations, 
and section 3c describes an Overdetermined Fixed Matrix (OFM) approach that is most 
desirable. Section 4 discusses tests of the OFM method using simulated sources with and without 
random errors in the AEs . Section 5 introduces a “Dimensional Reduction” (DR) method that 
uses the OFM approach to effectively reduce the number of unknowns in a general eight 
parameter (two-charge) model to just four unknown parameters. Section 6 uses the DR method 
to analyze over one hundred lightning flashes in a thunderstorm at the KSC-ER and then 
discusses the advantages and limitations of the DR method. 

As an alternative to the DR method, section 7 introduces a simple charge retrieval 
algorithm that uses KSC-ER Lightning Detection And Ranging (LDAR) data to constrain the 
charge locations, and a linear least-squares matrix inversion to determine the charge values. 
Section 8 applies this (computationally fast) method to derive charge solutions for a variety of 
Florida thunderstorms, and compares the results to independent satellite- and ground-based 
lightning observations. Finally, section 9 summarizes this work. 


2. Goodness-of-Fit Analyses 

A method for retrieving the optimum parameters of an assumed lightning charge model 
by numerically minimizing a nonlinear, chi-square error function was introduced by Jacobson 
and Krider (1976). Given m field change observations at the ground (i.e., A E n i = \,...,m) the 
reduced chi-square function is 


* 2 (a) 


1 ^ (M,.(a)-A£,) 2 

m~n t; erf 


( 1 ) 


4 



where M ; .( a) is the model field change at the i b site, and cr i is the observation error in AE i . The 
model M t is a function of n unknown model parameters that are contained in the n-vector 
a = (a p ..., a n ) . The difference (m - n) is termed the “degrees of freedom.” 

Simple examples of the model M,( a) , and the ones of greatest interest here, are the one- 
and two-charge models. Denoting the one-charge model by the subscript “1”, it can be written 




Z0£ 

InSoi^-xf+^-yf+^f 2 ’ 


( 2 ) 


where (X , , Y j ) is the location of the observation AE i , and a, =(x,y,z,Q) . This model assumes 
that the ground is a flat conducting plane of infinite extent. We use the standard physics 
convention so that M u is the change in the vertical electric field at the i' h site; hence, a 

monopolar cloud-to-ground lightning that transfers negative charge to ground would effectively 
deposit a single positive charge ( Q>0 ) at location ( x,y,z ) above the conducting plane, and 

thereby results in a negative change in field, M u < 0 . Using superposition and (2), the 8- 
parameter two-charge model can be written as 


M 2 ,(a 2 ) = 


-1 

2 


Q z 


Q'z 


[(W,. -x) 2 + (Y i -yf +z 2 f n [( X , -x'f +(Y, -y'Y +z' 2 ] 


f\2 , f2 n3/2 


( 3 ) 


where a 2 =(x,y,z,Q,x ,y ,z ,Q) . 

According to the statistical theory discussed in Bevington and Robinson (2003, Chapter 
1 1), a value of % 2 = 1 implies that the model-fitting function M, (a) in (1) is a good description 


of the observations provided that the values of cr j are accurately estimated. 
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To understand how well one must model a given set of A Es so that % 2 = 1 , suppose that 


the error in A E is constant, that is, cr. =<r; i = 1 . In this case (1) becomes 

i m 

Z 2 (a)= £(M,(a)-A E t f , 

(m-n)cr 


(4) 


and a “good fit” (i.e., % 2 = 1 ) implies that 


2 (Mi (a) - AE, ) 2 - (m-n)cr 2 . 

/=! 


(5) 


In addition, note that the rms modeling error, R (in V/m), is 


R = 


m 


£(M,.(a)-AE,.) 2 



( 6 ) 


Substituting the right hand side of (5) into (6) gives the approximate value of R required to 
obtain a “good fit,” that is, 

R = £cr => "good fit" ; 1 <n< m , (7) 


where E, = (1 -n/m) in . For the m = 31 sensors in the field measuring network at the KSC-ER, 
E, takes on the values shown in Table 1. Note that R must always be less than a to be viewed 
as a good fit; i.e., 0 < E, < 1 always holds. In addition, note that E, decreases as n increases. This 
means that a complicated model must describe the A Es more accurately than a simple model in 
order to be regarded as a “good fit.” For example, the one-charge model (« = 4) has E, = 0.933 , 
but the two-charge model ( n = 8) has E, = 0.861 . The improvement in modeling accuracy that is 
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required to maintain a “good fit” is not linear in n. As n increases, £ decreases at a progressively 
faster rate because the derivative d^jdn is inversely proportional to £ , i.e. 


i 



( 8 ) 


In summary, when applying models of increasing complexity, one must be sure that they 
can reduce the value of R by an amount that is sufficient to obtain a good fit. In the limit as 
n -> m , the model must describe the A Es perfectly in order to achieve a good fit. Again, rather 
than increasing the model complexity, n, Occam’s Razor suggests one should reduce model 
complexity as much as possible while still trying to maintain a good fit. 


3. Analytic Solution to the Point Charge Model 

a. Four station solution 

Appendix 1 of Krehbiel et al. (1979) describes an analytic method for determining the 
four unknown parameters of the Q-model from four ground-based observations of A E . This 
approach generalized an earlier method due to Fitzgerald (1957) that had certain limitations. By 
replacing the model field changes M Xi in (2) with the (absolute value of) the real field changes 

A E t produced by a monopolar flash, one obtains 



_2jf_ 

2 ^ 0 [(A. - a -) 2 +(Y i -y) 2 +z 2 ] 3/2 


( 9 ) 


Rearranging terms, and then differencing pairs of the i = 1,...,4 equations, Krehbiel et al. (1979) 
obtained 
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2(X i -X J )x + 2(Y i —Yj^y + UjjTj = X? +Y, 2 -X 2 -7/ , 


( 10 ) 


where U tj =[(4 ne o ) 2/j (| A E, | lh -| A£\ | 2/3 ), and rj = (2 \q\z) 2 ' 3 . [Note: Krehbiel et al. actually 
omit all the absolute value signs shown here, thereby restricting their result to the case Q > 0 . 
The results shown here are valid for any polarity of Q . ] Taking a combination of 4 sensors 2 at 
a time gives 4!/(2!2!) or 6 possible combinations. Krehbiel et al. pick 3 of these 6 to obtain 3 
equations in 3 unknowns, i.e. 

~2{X,-X 2 ) 2(7, —Y 2 ) U u Jxl [X?+Y*-Xl-Y*' 

2(X l -X 3 ) 2(7, -Y 3 U x 3 y = X 2 +Y 2 -X 2 -Y 2 , (11) 

_2(X,-X 4 ) 2(7, -7 4 ) [Z, 2 + 7, 2 -X 4 2 - 7 4 2 _ 

or in matrix-vector notation 

Av = b . (12) 

The column vector solution for the unknowns is v = col(x,y,rj) = A"'b . Substituting the (x,y) 
derived from this solution into (9) and noting that \Q\z -\rj 3 ' 2 , one can determine z . With z 
and the definition of 7 one obtains \Q\ . The sign of Q is the same as the sign of the A Es, hence, 

the problem is solved. Note that this method can only be used when the A Es have the same 
polarity, i.e., the solution is only possible for a monopole source. 

b. Overdetermined Variable Matrix (OVM) approach 

As discussed above, Krehbiel et al. (1979) used four ground-based A E observations to 
obtain the four unknown parameters ( x,y,z,Q ) in the Q-model. However, retrieval accuracy can 
be improved if all m observations are used at the same time. This assumes that the observation 
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errors are all similar in magnitude, and that the distances between the observation sites are 
sufficiently large so that the A matrix is not ill-conditioned. Inversion of an ill-conditioned 
matrix produces large retrieval errors because small eigenvalues magnify observation errors 
(Twomey, 1977; Chapter 6). 

Instead of extending their analytic solution to analyze all m observations of A E at the 
same time, Krehbiel et al. (1979) computed several four-station solutions using different 
combinations of sensors, and then averaged the results. This approach is unconventional in 
inversion theory (Twomey, 1977), and it is also computationally intensive. For m- 31 sensors 
(the number of field mills in the network at the KSC-ER), there are -4)!4!] = 3 1,465 

four-station solutions to compute. It is preferable to extend the difference scheme in (11) so that 
all m observations can be analyzed at the same time; this gives the overdetermined system of 
equations 

~2(X 1 -X 2 ) 2(7, —Y 2 ) U n lr -I I" Xt+Y?-Xl-Yf 

. . , y = ’ • ( 13 ) 
7iX x -X m ) 2(7, -Y m ) U h JH [ X]+Y?-Xl-Yl_ 

Defining again the matrix on the left hand side as A , and the vector on the right hand side as b, 
one can pre-multiply (13) by A (the transpose of A) and invert to obtain the solution 

v = (AA) _1 Ab . (14) 

This is a least-squares solution (Twomey, 1977; Chapter 6), and the optimum values of 
( x,y,z,Q ) can be computed from v as discussed in section 3a. Again, since the solution (14) 
uses all m observations to obtain the four model parameters, it is expected to provide a more 
realistic answer than any single solution derived from (12). Furthermore, external constraints can 
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be added to the solution process if desired (see Twomey [1977] for an introduction to 
constrained linear inversions). 

c. Over determined Fixed Matrix (OFM) approach 

The difference equations (10) were formed in order to remove the nonlinear terms x 2 ,y 2 , 

and z 2 ; however, these differences are not required to solve the problem analytically. An 
alternate and more practical solution (that we will ultimately use in section 5 to follow) can be 
found. We begin with (9), i.e., we assume that the source producing the A Es is a monopole so 
that all the field changes have the same sign. Raising (9) to the power gives 

d i = M{(X i -x) 2 +(Y i -y) 2 +z 2 ] 

= (X? + Y 2 )w- 2 X t wx - 2 Y t wy + wr 2 , 


where 


d, = A E, 


- 2/3 


W = 


2 7T£ c 

\Q\ Z 


\ 2/3 


2 2 2 2 
r = x +y+z 


For i = \,...,m sensors, (15) yields the (m x 4) system 


d, 

m 


'(X, 2 + if) -2X, 

-27, 

r 

w 




wx 




wy 

(X 2 m +Y„ 2 ,) ~2X m 

-2K, 

i 

wr 2 


that can be written in matrix-vector form as 

d = Ts . 


(16) 


(17) 


(18) 
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Premultiplying (18) by T and inverting gives 


s = Nd , (19) 

where N = (tT)“'T is termed the “network matrix” because its elements are determined only 
by the geometry of the sensor network. The model parameters can be obtained from the elements 
of s = (s l ,s 2 ,s 3 ,s i ) as follows 

x = s 2 /s l 

y = s 3 / s i 

z = (s l s 4 -st-s*y l2 /s l ^ 20) 

Q - (2^ 0 AE 1 )/[(^ 4 - S / 2 - Vs f ■ KQ . 

Note that the sign of Q is the same as any of the AEs, for example AE l /\AE l \ as shown in the 
last equation of (20). 

A distinct advantage of the solution in (19) over that provided in (14) is that the network 
matrix N need only be computed once (for a given network geometry), whereas the matrix A 
must be computed for each lightning flash because it contains the A E observations. Since the 
retrieval accuracy of a matrix inversion depends on the mathematical condition of the matrix, it 
is clearly better to have this condition fixed (and well understood), as in the case of N , rather 
than have it vary from flash to flash as with A . Of course, if any sensor fails or is removed from 
the analysis for any other reason, then N will change. 

In summary, (20) provides a way to use all m A E observations to infer the optimum 
source parameters ( x,y,z,Q ). Whereas the method discussed in section 3b obtains a solution by 

removing the nonlinear terms ( x 2 ,y 2 , and z 2 ) through differencing, and results in a variable A 
matrix from flash to flash, the OFM method provided in this section demonstrates that a solution 
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can be obtained without having to remove the nonlinear terms. The OFM approach eliminates 
the variable A matrix altogether, and effectively replaces it with the “fixed” network matrix N . 
This matrix only has to be computed once for any given network geometry, and it is not difficult 
or time consuming to compute because it involves the inversion of one small (4 x 4) matrix TT . 
After N is computed for a particular network, no other matrix inversions are required to 
retrieve monopolar flashes detected by that network. All the A E values are contained in the d 
vector and each monopolar solution is obtained from the productNd = s , followed by application 
of (20). 

4. Simulated Tests of the OFM Approach 

Before using the OFM method to help analyze natural lightning, we tested it on a number 
of simulated monopole sources that were placed over the sensor network at the KSC-ER. With 
no errors in the simulated AE observations, the OFM method retrieved the parameters of an 
arbitrary monopole exactly (apart from very small round-off errors). When a 5% random error 
was added to each of the simulated A E observations, the A E modeling errors and solution 
retrieval errors shown in Figure 1 were obtained. 

To cover the region in Figure 1, we placed sources on a horizontal grid that had 
dimensions 59 km x 59 km (resolution 1 km) and an altitude of 6.5 km. The error statistics were 
generated by analyzing 100 monopole sources (each 25 C) at each grid point. To obtain a set of 
simulated A E observations from a source, we computed the field at each site due to the source, 
and then added a random error to each field value. The random errors were chosen from a normal 
distribution that had a mean of zero, and a standard deviation equal to 5% of the field value. The 
percent AE modeling error, P, given in Figure la is defined as P =100xjAE-Mj/|AE| , where 
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AE and M are m-vectors of the simulated and model field changes, respectively. Using the key 
in Figure 2, note that horizontal and vertical errors in the retrieved distances are only 0.1 to 0.5 
km over the network, and the errors in Q are less than 1 C over a substantial portion of the grid. 
For the most part, the errors in horizontal distance, altitude, and Q increase outside the network. 
This behavior is typical for all inversion algorithms (see, for example, the tests of the multipole 
method in Appendix A of Koshak et al., 1999). In addition, note that it is more appropriate to 
plot P (Figure la) rather than the reduced chi-square. This is because the retrieval errors are 
directly calculated in these simulations, so the chi-square statistical theory for rating solution 
quality and estimating retrieval errors is not needed. Nonetheless, P and j 2 are mathematically 
dependant variables. For an assumed fixed a equal to 5% of the rms AE value for a monopolar 
flash, the green regions in Figure la (P ~ 5% to 10%) correspond to a range in ^ of 1.15 to 

4.59 which typically indicates physically reasonable fits. 

Because the OFM approach is essentially a single matrix inversion method for retrieving 
monopoles, it (reasonably) requires that all A Es be nonzero and of the same polarity. In addition, 
errors in AE can cause the difference (^ 4 -[s\ + s ] ]) in (20) to be negative, and this produces a 

complex number for the altitude, z. Errors in AE can also make the difference too small or too 
large, thereby producing a value of z that is erroneous based on physical grounds. Hence, if any 
of the above conditions (one or more zero A Es, bipolar AEs, a complex or unphysical z) occur, 
that test source has been ignored and another source has been generated as a replacement. To 
generate the statistics in Figure 1, we had to replace 96,179 sources. That is, we analyzed a total 
of 456,179 monopoles, removed 96,179, and this leave a total sample size of 360,000 (100 
monopoles times 3600 grid points) for the statistics in Figure 1. 
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In summary. Figure 1 shows that the OFM method provides a nice way to analyze the 
field changes produced by a single point (or spherically symmetric) change in the cloud charge 
distribution, but because A E observation errors (random and/or systematic) can “rupture” this 
pattern, the OFM method cannot invert all cases. However, in those cases where random errors 
are the primary cause of such a rupture, it is possible to add artificial random errors to the 
simulated observations in an attempt to remove some of the random errors that are responsible 
for the rupture. This removal process or “monopole unmasking” is an unorthodox technique 
because one is combating errors by adding more errors, however we have found it to be very 
useful. For example, we have re-run the same monopole simulation, but this time whenever we 
encountered a ruptured set of simulated A E values, we added a random 6% error to them, and 
then analyzed them again using the OFM method. If the pattern was still ruptured, we again 
added a 6% error to the (original) simulated measurements and re-analyzed them again with the 
OFM method. This process was repeated (up to 100 times) to see if the OFM method could 
eventually find a solution. We call this process the “extended OFM method,” and Figure 3 shows 
the results of this approach. The retrieval errors in Figure 3 are generally similar to those in 
Figure 1 (which indicates that adding errors can be useful in finding the hidden symmetry of 
monopole patterns). Whereas the standard OFM method could not analyze 96,179 (i.e., 
96,179/456,179 = 21.08%) monopoles due to ruptures, the extended OFM method had only 198 
(i.e., 198/360198 = 0.05%) that it could not analyze. 
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5. Dimensional Reduction (DR) Method 


a. Model complexity versus model capability 

One can define the complexity , n , of a given charge model to be the number of unknown 
parameters that it contains. For example, the one-charge model in (2) has n = 4 because it 
contains four unknown parameters ( x,y,z,Q ). By comparison, the two-charge model in (3) has 
n = 8 and is rather complex. As n increases, the model must do a progressively better job of 
describing the A Es in order to obtain a “good fit” (see (7) and Table 1). If a model has the same 
number of parameters as there are measurements; i.e., n = m , the model must describe the AEs 
exactly in order to achieve a “good fit.” So two disadvantages of a complex model are that it 
must meet stringent requirements to obtain a “good fit,” and it does not obey Occam’s Razor. 

The capability , k, of a particular model is the number of charge parameters that the 
model can describe. For example, the one- and two-charge models have k = 4 and a: = 8, 
respectively and for these two models, k = n . In general however, k need not always equal n . 
For example, the OFM method (see section 3c) can be used to reduce model complexity without 
reducing the model capability, and we term this approach “dimensional reduction” (DR) 
because the number of unknowns in the chi-square error function is reduced (see section 5b for 
more details). 

Given the above definitions, one desires to use a model that has minimum complexity, 
but sufficient capability. Because the one-charge model has just 4 unknown parameters, it cannot 
model many lightning events (e.g., intracloud discharges). Conversely, the two-charge model can 
describe many types of events (cloud-to-ground and intracloud discharges), but it has a high 
degree of complexity ( n = 8) . In between these cases, there is the 6-parameter 

(x,y,z,p xi p y ,p z ) point dipole model (Fitzgerald, 1957), the 7-parameter (x,y,z,x',y',z',Q) 
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discrete dipole model (Wilson, 1916), and the 7-parameter (x,y,z, Q,p x ,p y ,p z ) multipole 

model (Koshak et al., 1999). All of these models are more capable than the one-charge model but 
less complex than the two-charge model. 

b. The DR model 

It is possible to construct a model with complexity n = 4 , yet a capability k = 8 . 
Consider an idealized lightning flash that deposits two point charges a,=(x,y,z,0 and 
aj s (x',y',z',Q ') , and that is perfectly described by the two-charge model in (3). Using (2), the 
A Es at the ground due to this flash can be written as A£j. = M l/ (a 1 )+M ll .(aj). Therefore, given 
a, and the AC values, the point source aj must be chosen such that M u (aj ) = AE j - M l| . (a, ) , 
where the difference, AE l , is termed the “residual.” In other words, the two charges are 

interdependent, and one can write aj =(x',y,z',0) = (x'(a 1 ),y(a 1 ),z'(a,),0(a 1 )), or in vector 
function notation, aj =aj(a,). This interdependence motivates the definition of the 4-parameter 
“DR-model,” Z>.( a^ , given by 


A(».) 


,) + T,.(a|), a,ef 
i A a,fSF 


( 21 ) 


where T(a,) is a composite function defined by 


Y^sM^aKa,)): 


-Q' (a,)z'(a,) 


2ns 0 [{X, -x'(a,)) 2 +0 -y(a,)) 2 +(z'(a,)) 2 ] 


3/2 


( 22 ) 


Hence, D^a,) represents the superposition of fields produced from two point sources, yet it 
contains only four unknowns. It differs from the standard two-charge model 
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M 2j (a 2 ) = M 2j .(ai,aj) given in (3) because the relationship a', = a'j(a,) has been used to reduce 
model complexity; section 5c below discusses how the relationship aj =a' 1 (a 1 ) is determined. 
Since D j (a l ) is only a 4-parameter model, a good fit can be achieved wheni? = 0.933cr , and this 
is easier to achieve than the value R = 0.861(7 associated with the 8-parameter two-charge model 
M 2j (a v a[) (see Table 1). Finally, /? is a very large constant that represents a physically 

unrealistic field change value that D i ( a,) is driven to if a, strays outside a physically reasonable 
vector space V ; see section 5d for further discussion on how aj is constrained. 


c. Finding the optimum point sources 

In the DR method, the optimum source parameters a, = (x, y, z, Q) are found by 
minimizing the reduced chi-square function 


X 2 ( »,) = 


1 a,)-Ag,.) 2 

m -n~f erf 


(23) 


We perform a grid search of a slab that intersects the negative charge region (“N-region”; 
Williams, 1989) as shown in Figure 4. The grid search is followed by a Powell minimization 
(Press et al., 1992) that iteratively adjusts the parameters in aj until the minimum is found. If a, 

strays outside the hypervolume V given in (21), then the value of j 2 increases dramatically (i.e., 
X 1 °c 0 1 for a, <£V). At each step of the numerical minimization, the parameters (x\y',z',Q') 
are retrieved using the OFM method, with the components of d redefined as follows 

d, = dfaf = \AE t -M 1 ,.(a 1 )j _2/3 . (24) 
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That is, instead of d t depending on A £), it now depends on the residual (A/?,. -M u ) . In effect, 
the field change pattern is split up into two parts: one part that is described by the source 
( x,y,z,Q ), and the other part (or residual field change pattern) that is described by the source 
(x',y',z',Q') . If the residuals are of mixed polarity (or have any other “ruptures” as discussed in 
section 4), then a ( will not be regarded as acceptable and/or the errors in A E will be viewed as 
pathological. 

Note that the DR method can also be used to effectively reduce a 3-charge (12 parameter) 
model down to 8 unknown parameters. In this case, d t = d i (n 2 ) = \AE j - M v (a 2 )| , where M 2 f (a 2 ) 

is the two-charge model in (3). In general, the DR method can convert any 7V-charge (4N 
parameter) model to a (47V-4 parameter) model by using the following definition 


dj — d i (a ( ^ N _ i)) — 


N - 1 




Qj z j 


t2ns 0 [{X- Xj y HY-y.Y 


(25) 


where a^_^ — x N _ x , y N _ { , z N _ x , Q N _ X ) . 

d. Constraint details 

Recall that we impose constraints on the source parameters (x,y,z,Q) = a 15 where 
a, e V , and V is an “acceptable vector space.” Our constraints are straightforward, but highly 
advantageous: the source location ( x,y,z ) is restricted to the N-region of a thundercloud, and 
the value of Q is restricted to values commonly observed or inferred for the N-region. 

Since most lightning flashes begin within or near the N-region of the thundercloud, there 
is little loss in generality by imposing these constraints. In fact, because the N-region usually has 
a small vertical extent and is less diffuse than the upper positive or P-region, and because Q > 0 
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normally holds in the N-region, these constraints substantially reduce the size of V . A smaller 
V means it is easier to minimize the chi-square function, and this translates into reduced 
computer time, and improved accuracy. Because it is difficult to accurately retrieve the charge 
parameters when the source is not located over or near the measuring network (Koshak and 
Krider, 1999), a horizontal range constraint can also be imposed. Such a constraint further 
reduces the size of V . 

Applying these ideas, and keeping in mind the results of prior investigations (e.g., 
Koshak and Krider, 1989), the following constraints (on both aj and aj) can be imposed when 

analyzing warm-season storms at the KSC-ER: 

(a) -30 km < (x and y)< 30 km 

(b) 5 km < z < 8 km 

(c) O<0<1OOC 

(d) -30 km < (x' and y') < 30 km 

(e) 4 km < z' < 16 km 

(f) 0 < \Q'\ < 100 C 

Again, if during the numerical minimization of the chi-square function, any element of aj 
violates any of the above constraints, where a', =aj(a,), the chi-square is assigned a large value 

proportional to 0 1 . The chi-square must also be assigned a large value if the residual field change 
pattern is ruptured (see section 4). For example, if the field changes are of mixed polarity, OFM 
retrievals of aj will not be possible, and the chi-square must be assigned a large value. 

Note that the restrictions on the altitude, z, in (b) above are fairly tight, but z can vary 
over a much larger range, (e). This is by design because the OFM method can cover a larger 
altitude range more quickly and more thoroughly than a numerical minimization procedure. For 
example, it would not be wise to initialize the a, source by placing it in the P-region of a 
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thunderstorm because we know that this region is likely to be larger than the N-region, and 
therefore it would take more numerical iterations to minimize the chi-square function, all else 
being equal. 

e. Solution refinement and associated errors 

After minimizing the chi-square, a (preliminary) two-charge solution given by 
a,c=( x c>yc’ z c>Qc) = ( r c>Qc ) and *4 = (x' c ,y' c ,z' c ,Q' c ) = (r' c> Q' c ) is obtained. To improve this 
solution, a cubical box (dimension 1 km) is centered on each of the charge locations (r c ,r c ') . We 
then sample 5000 locations within each box; one sample is the center of the box, and the 
remaining 4999 samples are randomly selected. From each pair of sampled locations 
(r ; ,r y '), / = !,..., 5000, we construct the kernel matrix K } =K(r\,r'), and compute the standard 
linear least-squares charge retrieval 

= (K.K y r‘K.AE; y = l,...,5000 ( 26 ) 

[see Koshak et al. (1999) for a description of the kernel matrix and linear charge retrieval 
theory]. Of the 5000 solutions { (Q, , r y , Q] , rj ) , j = 1,...,5000 }, the solution that produces a 

minimum in the chi-square is saved as the optimum solution. Note that =r c and vl = r ' c ; i.e., 

the value j = 1 corresponds to the center of each box. However, Q c and Q' are from the chi- 
square minimization process using the OFM approach and therefore are not equivalent, in 
general, to the respective linear least-squares values Q x and Q[ , 

Because a linear least-squares solution (26) has been used to find the charge values, 
standard linear least-squares estimation theory can be used to estimate the errors in these values. 
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Specifically, the covariance matrix C = (KT 'K) 1 (also known as the inverse curvature matrix) 
contains diagonal elements that are the variances of each of the retrieved charges, and the square 
root of these variances provide estimates of the charge errors. Here, T is a (m x m) diagonal 
matrix whose diagonal elements are the squared observation errors a] , i = 1 . Of course, 

error estimates derived from C are only valid if the charge locations themselves are reasonably 
accurate (and the model is perfect). 

6. Application of the DR Method 

We have used the DR method to analyze a 30 minute interval (20:30-21:00 UTC) of a 
warm-season Florida thunderstorm that occurred at the KSC-ER on 21 September 1998. Sensors 
#1 and #32 were faulty during this period and were not included in the analysis. Flashes that 
produced a A E> 100 V/m at 2 or more sites were analyzed; i.e., a total of 154 flashes. Of these, 
115 flashes produced a reduced chi-square less than 5; the average reduced chi-square for these 
115 solutions was 1.058 with a standard deviation of 1.017. [We used a chi-square threshold of 5 
because there is uncertainty in the value of a in (4), and because many solutions meeting this 
criterion were physically reasonable.] 

The altitudes of the 115 charge solutions (large circles) are plotted in Figure 5 together 
with LDAR data (small dots). Positive charge deposition is denoted with red circles, and 
negative deposition with blue. Larger circles indicate greater amounts of charge deposition; i.e., 
the radius r of each circle has been scaled to the vertical axis (see Koshak and Krider, 1 989) and 
is a solution to the breakdown field equation E b = Qj (4ns 0 r 2 ) with E b = 400 kV/m. The LDAR 

system maps the source locations of the VHF radio impulses produced by both intracloud (IC) 
and cloud-to-ground (CG) lightning with a time-resolution of about 100 ps and a spatial accuracy 
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of hundreds of meters or less for flashes that are within or near the perimeter of the KSC-ER 
field mill network (Maier et al., 1995). For lightning that is within about 40 km of the KSC-ER, 
the LDAR location accuracy is better than 1 km (see Boccippio et al., 2001, Parts 1 and 2 for 
additional details on LDAR performance). The LDAR data are color-coded according to time in 
the plot projection. 

To understand the quality of the charge solutions in Figure 5, we have examined 
individual plots of each flash to see how well the lightning charge locations and the LDAR data 
agree. Given that the A E and LDAR data are completely independent datasets, it was 
encouraging to find that a high percentage of the model solutions had (x, y) locations that were 
within 1 km of, or co-aligned with, the ( x,y ) positions of LDAR sources. Overall, we found that 
there were three types of solutions: (I) those that lined up well with LDAR sources in x, y, and z, 
(II) those that lined up well in x and y but had offsets in z , and (III) those that had offsets in x, 
y, and z. 

Figures 6, 7, and 8 show examples of solution types I, II and III, respectively. Most 
offsets in z (Figures 7 and 8) were caused by a fundamental “dipole ambiguity.” For example, 
consider a vertically oriented cloud discharge that produces a dipole moment change (20C)(3km) 
= 60 C-km. If the charge transfer is reduced to 15C and the discharge length is increased to 4 km, 
the dipole moment does not change; i.e., (15C)(4km) = 60 Ckm, and it is possible (depending on 
source location and network geometry) that the two A E patterns associated with these two 
sources are nearly indistinguishable. Moreover, any small features that distinguish the patterns 
can be completely removed by errors in A E . Hence, fundamental ambiguities arise whenever 
different dipole sources effectively produce the same set of A E values. For all the solutions in 
Figures 6-8, the charge retrieval errors (i.e. the square root of the diagonal elements of the 
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covariance matrix C) were typically less than 1 C. Of course, this magnitude of error applies 
primarily to type I (Figure 6) where the charge locations are reasonably accurate and the linear 
least squares error estimation theory can be used. 

Of the 115 flashes, 21 (18.3%) were exceptional fits (type I), 61 (53.0 %) could be fit 
reasonably well (type II), and 25 (21.7%) were type III. Additionally, 5 (4.3%) were complicated 
flashes that were too difficult to categorize into any one of the three main types, and 3 (2.6%) 
were false triggers in the electric field dataset (i.e. the A Es were not associated with LDAR 
events). Since in this case study a total of 82 (71.3%) could be fit either reasonably or 
exceptionally well, we believe that the DR method is a valuable tool for describing lightning 
field changes. 

In section 7 to follow, we will describe how we use LDAR data to constrain the charge 
location model parameters so that, by definition, they do not contain offsets in position. These 
charge retrievals will then be compared to both satellite- and ground-based lightning 
observations in section 8. 

7. LDAR-Constrained (LC) Method 

Because lightning effectively deposits charge within thundercloud charge centers and 
because LDAR traces the geometrical development of the lightning channel with high precision, 
the LDAR data provides an ideal constraint for finding the best model charge solutions. In 
particular, LDAR data can be used to help determine both the horizontal and vertical positions of 
the model charges, thereby eliminating dipole ambiguities. 

The LDAR-constrained (LC) charge retrieval method is straightforward. LDAR data for 
individual flashes are displayed on the computer monitor, and a mouse is used to click the 
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expected locations of the charge centers (usually the endpoints of the LDAR-mapped lightning 
channel). For example, a typical vertical or slanted cloud discharge between the N- and P- 
regions would receive two mouse clicks, one at the lower end of the channel within the N-region, 
and one at the upper end within the P-region. Calling these endpoint locations r, and r,' , 
respectively, we center a cubical box (dimension 1 km) about each of these locations and 
internally sample the box volumes exactly as described in section 5e. That is, we compute 5000 
solutions {(Qj,rj,Qj,r'j) , j = 1,...,5000}, and the solution that produces a minimum in the chi- 

square is saved as the optimum solution. Hence, the discharge orientation is obviously tightly 
constrained by the LDAR channel endpoints, but the model charges are allowed to deviate a bit 
from these endpoints in order to obtain a better fit to the A E observations. 

8. Satellite-Observed Flashes 

The results of the LC method have been compared to the locations of optical pulses 
detected by the Lightning Imaging Sensor (LJS) on the low Earth-orbiting Tropical Rainfall 
Measuring Mission (TRMM) satellite (Christian et al., 1999). Several flashes from different 
storms that occurred near the KSC-ER have been studied in detail. To augment this investigation, 
we have also compared both datasets to data provided by the U.S. National Lightning Detection 
Network™ (NLDN) and the Cloud-to-Ground (CG) Lightning Surveillance System (CGLSS) that 
covers the region of the KSC-ER. Of course, the charge solutions have also been compared to 
the LDAR data that have been used to constrain the charge locations. The CGLSS and NLDN 
systems locate the ground strike points of the individual return strokes in CG flashes with an 
accuracy of a few hundred meters (see, for example, Cummins et al., 1998a,b). Also, Thomas et 
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al. (2000) have found good overall agreement between Lightning Mapping Array (LMA, a 
system similar to LDAR) sources and LIS flash locations in Oklahoma. 

Although it is rather rare to have a TRMM overpass at the same time a thunderstorm is in 
progress at the KSC-ER, 7 such periods have been identified within the years 1998-1999 and 
were analyzed for this study (see Table 2). Each overpass lasted from 2 to 3 minutes, and the 
total observation time was 1 5 minutes. The KSC-ER field mill (FM) network detected a total of 
94 discharges (each producing AC > 100 V/m at 2 or more FM sites) during this 15 minute 
period, and the LDAR system detected all 94 of these discharges. Figure 9 shows examples of 
two overpasses from the unique perspective of the LIS low-Earth orbit; both the LIS background 
image in the near infrared and the density of LIS-detected flashes are provided. 

The CGLSS data were not available for days 180, 219, and 229, but when they were 
available, there were always more CG events in the CGLSS dataset than were detected by the 
NLDN. Also, many of the 94 discharges were cloud flashes that were not reported by the CGLSS 
or NLDN unless they were incorrectly classified. 

The lightning flashes that were not reported by LIS were divided into two categories, a 
Legitimate Miss (LM) or a Not Processed (NP) type. Typically, the LMs were due to LIS data 
post-filtering, “pixel splitting” of radiance, and/or cloud attenuation. The NPs were due to LIS 
buffer overflows, flashes (just) outside LIS field-of-view, and spacecraft incidents (e.g. telemetry 
bit errors, and spacecraft attitude maneuvers). Of the 94 discharges in our case study, 93 appear 
to be regular “flashes,” and one was a small breakdown event that preceded a normal CG flash. 
Of the 93 flashes, LIS did not report 17, and of these, we estimate that 6 were LMs, and 1 1 were 
NPs. The NPs were mostly due to buffer overflows and/or outside field-of-view cases, although 
one flash was probably lost due to a bit error in the spacecraft telemetry. 
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Lightning charges were computed for all 94 discharges and compared to the patterns of 
AE obtained from the field mill network, the spatial-temporal development of the LDAR 
sources, with the locations of LIS optical events, and with the NLDN and CGLSS ground strike 
points, when available. As expected, the LDAR-constrained charge retrievals had locations that 
were in very good agreement with data obtained from the other lightning detection systems. With 
an assumed <j equal to 15% of the typical (rms) field change value, the median and average 
reduced chi-square for the 94 flashes was 3.69 and 9.54, respectively. A total of 55 (~ 60 %) of 
the 94 flashes had reduced chi-square values less than 5, and were good solutions. We also found 
physically reasonable solutions for larger reduced chi-square values (see Figures 11 and 13 to 
follow). Additionally, it is important to note that when one adds external constraints (e.g. LDAR 
data) to the retrieval process, it is common to have an increase in the chi-square value. This is 
because the (now constrained) model parameters can no longer be adjusted in any way desired to 
mathematically minimize the chi-square. Therefore, when applying valid external constraints 
(such as LDAR data), it is incorrect to infer that any increases in chi-square imply a lower 
quality solution and/or a physically unrealistic model, especially when AE errors are actually 
large when assumed to be reasonably small. 

Figures 10-15 show detailed results for 6 of the 94 discharges that were analyzed using 
the LC method. The red and blue circles show the locations of the positive and negative 
lightning-charges, Q and Q ' , respectively, that the lightning effectively deposited in the 
thundercloud. The values of Q and Q' are given at the top of each plot in Coulombs. The 

curved lines in the plan views on the lower left show the contours of AE in V/m; red contours 
indicate positive values, and blue contours indicate negative values. The reduced chi-square 
values associated with the solutions in Figures 10-15 are, respectively: 1.82, 11.8, 0.93, 6.55, 
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0.79, and 4.4. The panels on the upper right show plots of altitude versus time (in GMT) with 
color coding; the total time is 2 seconds in each plot except for Figure 1 1 which is 1 second. The 
large squares show the ( x,y ) locations of the LIS optical events, and the large diamonds are the 
LIS flash locations. Since the LIS sensor does not directly observe the altitude of a flash, the 
optical events and locations are placed near the top of each altitude cross section plot for 
optimum readability, and to remind the reader that these data are derived from a satellite that is 
viewing the KSC-ER from above. The small dots show the locations of LDAR sources, and each 
dot has been color-coded with respect to time using the same color sequence as in Figure 5. The 
small squares show the locations of the first LDAR source (black) and last (orange) in the 
discharge. The black and orange “X” symbols show the locations of negative ground strike 

points from NLDN and CGLSS, respectively; a “+” symbol indicates that the ground flash had a 
positive polarity. 

The discharges in Figures 10-15 were selected because they include both CG and cloud 
flashes, and because they occur to the north (Figure 10), south (Figure 11), east (Figure 12), west 
(Figure 13), and over (Figure 14) the network. The discharge in Figure 14 is particularly 
interesting because it is a small breakdown event that precedes a CG flash. Even though the flash 
types and locations are diverse, the LDAR dataset provides good location constraints so that the 
linear least-square charge retrievals are physically reasonable in both magnitude and polarity. 

We believe that the discharge in Figure 1 1 is a cloud flash, rather than a positive CG as 
indicated by NLDN/CGLSS. This belief is supported by several factors: (1) the charge retrievals 

t sjc 

are roughly equal and opposite; i.e., a 10% charge retrieval error is possible , (2) no LDAR 

* Taking the charge values literally implies that almost 4 C (40.070 C - 36.123 C) of 
negative charge is deposited to ground, thereby contradicting the polarity inferred by 
NLDN/CGLSS. 
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sources are below about 5 km, and (3) there are similarities in the field change characteristics of 
positive CGs and normal vertical cloud flashes. Previous investigators have recommended 
classifying positive flashes with peak currents less than 10 kA as cloud discharges [Cummins et 
al., 1998b; Zajac and Rutledge, 2001]; and a 15 kA threshold has been more recently used 
(Cummins personal communication, 2006). Since the peak currents inferred by NLDN and 
CGLSS are 15.3, and 14.0 kA, respectively, for the flash in Figure 1 1, this is one example where 
application of a simple 10 kA threshold criterion would appear inadequate. 

An analysis of the LIS raw data shows that the LIS instrument likely detected the small 
discharge in Figure 14, but that it was later removed by the LIS data processing algorithm. 
Because this discharge was less than 0.9 s before, and within 10 km of, the CG flash provided in 
Figure 15, we think that these discharges might be physically correlated to one another. From 
this point of view, the small discharge in Figure 14 is not considered a “legitimate flash” on its 
own, and one can argue that LIS properly reported just the primary portion of the overall 
phenomenon. 

9. Summary 

We have introduced an Overdetermined Fixed Matrix (OFM) method for finding the 
optimum parameters of a 4-parameter Q-model. This method offers two distinct advantages over 
a previous 4-station solution scheme: first, one can use more than four A E observations to 
obtain a highly constrained source retrieval. Second, the inverse matrix is “fixed” which means it 
does not depend on the values of A E (as in the previous 4-station solution scheme), and hence 
the matrix does not change from flash to flash. Because the matrix only depends on the 
measurement network geometry, only one matrix inversion is required for any given network. 
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This speeds up computations, and it also allows one to more thoroughly understand the retrieval 
errors associated with any particular network geometry since the eigenvalues of a fixed matrix 
are fixed (see Chapter 6 of Twomey, 1977 regarding the importance of eigenanalyses in linear 
retrievals). 

We have shown that the OFM method accurately retrieves the parameters of simulated, 
single-charge sources in the presence of A E errors. For 25 C monopoles placed 6.5 km over the 
field mill network at the KSC-ER, the horizontal ( x,y ) and vertical (z) retrieval errors were each 
less than 0.5 km and the charge retrieval error was less than 1 C. As with most analysis methods, 
retrieval errors increase for sources located off the network, and continue to increase with source 
range. 

Most importantly, the OFM method has been used to reduce the complexity (number of 
model parameters) of a two-charge model. This “dimensional reduction” (DR) approach 
significantly reduces the computational burden because the number of unknown parameters in 
the chi-square function is reduced from 8 to just 4. Because of this reduction in dimension, there 
is a reduced chance of getting “stuck” in a local minimum of the chi-square function (because 
there are fewer unknowns), and consequently the DR method has the potential to significantly 
reduce retrieval errors. 

One should note that the DR method is not simply a practical approach for solving charge 
retrieval problems. Rather, it introduces a fundamental, yet subtle, aspect of lightning charge 
modeling: the distinction between “model complexity” and “model capability.” Model 
complexity is the number of unknown charge parameters contained in the model, whereas model 
capability is the number of charge parameters the model can describe. [For example, the model 
provided in (21) has 4 unknown charge parameters (complexity = 4), but it can describe 8 charge 
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parameters (capability = 8).] We also showed that when the complexity of a model increases, one 
is forced to model the A Es more accurately in order to obtain a good fit (see (7) and Table 1). 
Since the DR method reduces model complexity but maintains model capability, it makes it 
easier to obtain a good fit. As an added benefit, the DR method also satisfies the simplicity 
principle of Occam’s Razor. 

Finally, we have used the DR method to analyze over one-hundred lightning flashes in a 
Florida thunderstorm, and it produced a high percentage of physically reasonable solutions, as 
judged by comparisons with LDAR data. Unfortunately, several DR solutions suffered from a 
fundamental “dipole ambiguity” that produces vertical offsets in the retrieved charge altitudes 
relative to those that are inferred from the LDAR dataset. This problem has been solved by 
using LDAR data to constrain the location of the lightning charges, and the resulting “LDAR- 
constrained” charge solutions have been compared to LIS optical sources and the locations of the 
ground strike points obtained from the NLDN and the CGLSS at the KSC-ER. We have given 
examples of these comparisons for several Florida thunderstorms, and in the future we expect 
that additional comparisons will provide relationships between lightning charges and the optical 
radiance, the number of LDAR sources, and in the case of CG flashes, the peak current and 
multiplicity. 
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Figure Caption List 


Fig. 1. Monte Carlo error analyses of monopole sources retrieved with the OFM approach 
showing (a) average percent A E modeling error, (b) rms horizontal distance error, (c) rms 
altitude errors, and (d) rms charge errors. The color key corresponding to these figures is shown 
in Figure 2. The error statistics are derived for 100 randomized sources at each grid point, with a 
grid resolution of 1 km. 

Fig. 2. The key applied to Figure 1 (and Figure 3 to follow) that associates color to error values 
for the three error categories shown. 

Fig. 3. Same type of error analysis as in Figure 1 except that a 6% random error was purposely 
added to the simulated AEs (in several independent trials) in order to nullify some of the 
simulated observation errors. This resulted in significantly increasing the number of monopole 
sources successfully modeled. 

Fig. 4. A summary of Dimensional Reduction (DR) where eight unknowns 
(x,y,z,Q,x',y',z',Q') are converted to just four ( x,y,z,Q ) using the OFM method. An example 
of one of the residuals is shown in purple. As the chi square function (23) is minimized, the 
source, a, (which produces field M u ( a,) at the ground), approaches both the location and 
magnitude of the lower thundercloud charge, and the OFM-retrieved source, a| , approaches the 
location and magnitude of the upper thundercloud charge. 
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Fig. 5 . Altitude vs. time plot of the DR-retrieved charges and associated LDAR data for a 30 
minute period of a thunderstorm that occurred on 21 September 1998. The red circles indicate 
positive changes in cloud charge, and the blue circles indicate negative changes. The small dots 
are the LDAR radio pulses (color coded with time). See Figures 6-8 below for further details on 
individual flashes. 

Fig. 6. Examples of charge solutions derived from the DR method (and picked from the results in 
Figure 5) that compare well with LDAR data. 

Fig. 7. Examples of charge solutions derived from the DR method (and picked from the results in 
Figure 5) that show charge location offset errors in the vertical relative to LDAR data. 

Fig. 8. Examples of flashes derived from the DR method (and picked from the results in Figure 
5) that show charge location offset errors in the horizontal and vertical relative to LDAR data. 

Fig. 9. Examples of two LIS overpasses of KSC on (a) 14 May 1999, and (b) 21 September 
1998. The background cloud image is not always collected if lightning rates are high. A total of 7 
overpasses (on seven different storm days) were analyzed in this study. 

Fig. 10 . A CG flash located north of the KSC-ER FM network. The red and blue circles denote 
lightning-caused changes in thundercloud charge given by Q and Q ' , respectively as derived by 
using LDAR constraints and linear least squares theory. The values of Q and Q' in coulombs 
are given at the top of the plot. The large squares are individual LIS optical events and the large 
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diamond is the LIS-derived flash location. The small dots indicate VHF source locations for the 


flash derived from the LDAR network; the small squares identify what VHF source occurred 
first (black) and last (orange) within the flash. The black and orange “Xs” mark the NLDN- and 
CGLSS- derived CG location, respectively. 

Fig. 11. An IC flash located well south of the KSC-ER network. See caption of Figure 10 for a 
description of the plot symbols. We believe this flash was incorrectly categorized as a positive 
CG by the NLDN and CGLSS systems (see main text for a discussion). 

Fig. 12. A CG flash (with IC component) located east of the KSC-ER network. See caption of 
Figure 10 for a description of the plot symbols. 

Fig. 13. An IC flash located slightly west of the KSC-ER network. See caption of Figure 10 for a 
description of the plot symbols. 

Fig. 14. A small discharge over the KSC-ER network that the LIS instrument detected, but that 
the LIS data processing algorithm removed. It precedes the CG shown in Figure 15 by less than 
0.9 s, and is also within about 10 km of the CG (see main text for more comments). The CGLSS 
and NLDN ground strike locations shown are associated with the CG that follows, but appear in 
this plot due to sensor timing errors. See caption of Figure 10 for a description of the plot 
symbols. 
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Fig. 15. The CG that quickly follows, and is likely correlated with, the small discharge shown in 


Figure 14. The CGLSS and NLDN ground strike locations for this CG appear in Figure 14 due to 
sensor timing errors. See caption of Figure 10 for a description of the plot symbols. 
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Values of the Coefficient 


Table 1. The numerical value of the coefficient £, in equation (7) as a function of model 


complexity n , when m = 31 sensors. The bolded values indicate the customary 1- and 2- charge 


models. 


n 

£ («; m = 

i 

0.984 

2 

0.967 

3 

0.950 

4 (e.g., the 1-charge model) 

0.933 

5 

0.916 

6 

0.898 

7 

0.880 

8 (e.g., the 2-charge model) 

0.861 

9 

0.842 

10 

0.823 

11 

0.803 

12 

0.783 

13 

0.762 

14 

0.741 

15 

0.718 

16 

0.696 

17 

0.672 

18 

0.648 

19 

0.622 

20 

0.596 

21 

0.568 

22 

0.539 

23 

0.508 

24 

0.475 

25 

0.440 

26 

0.402 

27 

0.359 

28 

0.311 

29 

0.254 

30 

0.180 

31 

0.000 
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Summary of Storm Periods 


Table 2. Summary of the Florida storm periods analyzed that were associated with LIS 
overpasses of the KSC-ER. 


Date 

Julian Day 

Time (GMT) 

# flashes 
(FM-derived) 

# Bad 
FMs 

21 Sep 1998 

264 

20:39-20:42 

13 

2 

08 May 1999 

128 

22:04-22:06 

11 

1 

14 May 1999 

134 

19:38-19:40 

22 

0 

11 Jun 1999 

162 

05:06-05:08 

14 

3 

29 Jun 1999 

180 

19:01-19:03 

21 

1 

07 Aug 1999 

219 

23:37-23:39 

7 

0 

17 Aug 1999 

229 

17:58-18:00 

6 

Total = 94 

0 
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X (km) X (km) 


Fig. 1 . Monte Carlo error analyses of monopole sources retrieved with the OFM approach 
showing (a) average percent A E modeling error, (b) rms horizontal distance error, (c) rms 
altitude errors, and (d) rms charge errors. The color key corresponding to these figures is shown 
in Figure 2. The error statistics are derived for 100 randomized sources at each grid point, with a 
grid resolution of 1 km. 
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30-50% 


>5km 

>2QC 

>50% 


Fig. 2. The key applied to Figure 1 (and Figure 3 to follow) that associates color to error values 
for the three error categories shown. 
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Fig. 3. Same type of error analysis as in Figure 1 except that a 6% random error was purposely 
added to the simulated A Es (in several independent trials) in order to nullify some of the 
simulated observation errors. This resulted in significantly increasing the number of monopole 
sources successfully modeled. 
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Fig. 4. A summary of Dimensional Reduction (DR) where eight unknowns 
(x,y,z,Q,x',y',z',Q') are converted to just four ( x,y,z,Q ) using the OFM method. An example 
of one of the. residuals is shown in purple. As the chi square function (23) is minimized, the 
source, a[ (which produces field M u ( a,) at the ground), approaches both the location and 

magnitude of the lower thundercloud charge, and the OFM-retrieved source, a \ , approaches the 
location and magnitude of the upper thundercloud charge. 
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time (GMT) 

Fig. 5. Altitude vs. time plot of the DR-retrieved charges and associated LDAR data for a 30 
minute period of a thunderstorm that occurred on 21 September 1998. The red circles indicate 
positive changes in cloud charge, and the blue circles indicate negative changes. The small dots 
are the LDAR radio pulses (color coded with time). See Figures 6-8 below for further details on 
individual flashes. 
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Fig. 6. Examples of charge solutions derived from the DR method (and picked from the results in 
Figure 5) that compare well with LDAR data. 
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Fig. 7. Examples of charge solutions derived from the DR method (and picked from the results 
in Figure 5) that show charge location offset errors in the vertical relative to LDAR data. 
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Fig. 8. Examples of flashes derived from the DR method (and picked from the results in Figure 
5) that show charge location offset errors in the horizontal and vertical relative to LDAR data. 
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Fig. 9. Examples of two LIS overpasses of KSC on (a) 14 May 1999, and (b) 21 September 
1998. The background cloud image is not always collected if lightning rates are high. A total of 7 
overpasses (on seven different storm days) were analyzed in this study. 
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May 14, 1999 (Day 134) 

Q=56.168, Q’ — 8.209 



Fig. 10. A CG flash located north of the KSC-ER FM network. The red and blue circles denote 
lightning-caused changes in thundercloud charge given by Q and Q ' , respectively as derived by 
using LDAR constraints and linear least squares theory. The values of Q and Q' in coulombs 
are given at the top of the plot. The large squares are individual LIS optical events and the large 
diamond is the LIS-derived flash location. The small dots indicate VHF source locations for the 
flash derived from the LDAR network; the small squares identify what VHF source occurred 
first (black) and last (orange) within the flash. The black and orange “Xs” mark the NLDN- and 
CGLSS- derived CG location, respectively. 
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June 11, 1999 (Day 162) 

Q=40.072, Q’ — 36.123 



Fig. 11. An IC flash located well south of the KSC-ER network. See caption of Figure 10 for a 
description of the plot symbols. We believe this flash was incorrectly categorized as a positive 
CG by the NLDN and CGLSS systems (see main text for a discussion). 
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June 29, 1999 (Day 180) 

Q=33.633, Q’=- 3.760 



x 








September 21, 1998 (Day 264) 

Q=2.4(52, Q’=-2.412 



Fig. 14. A small discharge over the KSC-ER network that the LIS instrument detected, but that 
the LIS data processing algorithm removed. It precedes the CG shown in Figure 15 by less than 
0.9 s, and is also within about 10 km of the CG (see main text for more comments). The CGLSS 
and NLDN ground strike locations shown are associated with the CG that follows, but appear in 
this plot due to sensor timing errors. See caption of Figure 10 for a description of the plot 
symbols. 
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